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SUMMARY 

A numerical method is presented which is valid for integration of the parabolic - 
elliptic Navier -Stokes equations. The solution procedure is applied to the three- 
dimensional supersonic flow of a jet issuing into a supersonic free stream. Difficulties 
associated with the imposition of free -stream boundary conditions are noted, and a coor- 
dinate transformation, which maps the "point at infinity" onto a finite value, is introduced 
to alleviate these difficulties. 

Results are presented for calculations of a square jet and var 3 dng -aspect-ratio rec- 
tangular jets. The solution behavior varies from axisymmetry for the square jet to 
nearly two-dimensional for the high -aspect-ratio rectangle, although the computation 
always calculates the flow as though it were truly three-dimensional. 

INTRODUCTION 

The calculation of free -mixing flows has> in the past, been accomplished through 
use of the boundary-layer assumptions in the two-dimensional or axisymmetric Navier- 
Stokes equations. The accuracy and validity of„these procedures rhaye been weir docu- 
mented in the literature (refs. 1 and 2). However, there are numerous situations where 
the flow cannot be considered either two-dimensional or axisymmetric. Jets issuing 
from rectangular orifices (see fig. 1), wakes behind any but the simplest bodies, and the 
flow downstream of a wingtip are examples of three-dimensional free -mixing flows where 
boundary -layer assumptions are invalid. The characteristic feature of these flows is the 
importance of diffusion in two spatial coordinates. 

These flows have certain characteristics which are common to boundary-layer 
flows; e.g., the velocities in the planes normal to the main-stream direction are usually 
much smaller than the main -flow velocity. Consequently, one expects gradients that 
exist in the cross directions to be larger than the gradients in the main -flow direction. 
Also, these are usually constant -pres sure flows. Therefore, it might be reasonable to 


use a boundary -layer scaling on the Navier -Stokes equations to effect some simplifica- 
tion; however, this yields an inconsistent set of equations. To lowest order, the cross- 
stream momentum equations reduce to statements that the pressure is coimtant, and the 
resulting number of equations is not sufficient to determine the remaining unknown 
quantities. , 

This paper presents the results of a method which overcomes this difficulty and 
allows numerical solutions of the parabolized Navier -Stokes equations. The method is . 
applied to a three-dimensional, s\q)ersonic, rectangular jet problem in which the ai^ect 
ratio is varied from one (square jet) to large values representative of slits. The range 
of applicability of the procedure is demonstrated from the near axisymmetry of the 
square, through a true three-dimensional region of moderate a^ect ratios, to a quasi- 
two -dimensional flow in the h^h-aspect-ratio limit which approximates a two- 
dimensional jet. 

SYMBOLS 


A,B 

Ca1[b].[c] 

D 

Dt 

h 

M 

Npr 

P 

R 

S 


transformation constants 

matrices defined by equation (15) 
vector defined by equation (15) 

convective derivative 

half -width of unit jet (reference length) 
Mach number 
Prandtl number 

pressure 
Reynolds number 
= SVTr 


S* Sutherland's constant 
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T temperature 

U velocity 

u,v,w component velocities in x-, y-, and z-directions, respectively 
x,y,z coordinate directions 

Ax,Ay,Az grid spacings in x-, y-, and z-directions, rei^ectively 
W vector of dependent variables (see eq. (16)) 

y ratio of specific heats 

A dilatation 

5 central difference operator 

transformed y- and z -coordinates 
viscosity 

p density 

i dissipation 

Subscripts: 

j,k discretized y- and z -positions, respectively 

r reference quantity 

* free stream 

- iq>per and lower grid spacings used with ponuniform grid 

Superscripts: 


I 


discretized x-position 



intermediate step of alternating -direction implicit method 
GOVERNING EQUATIONS 
Derivation of Equations 

The full three-dimensional Navier-Stokes equations are elliptic in character. The 
core storage available on present computers is insufficient to practicably handle any but 
the coarsest computational grids. Thus, methods to reduce the equations to a form more 
tractable for computation must be employed.- A true boundary -layer scaling cannot be 
used since it yields an inconsistent set of equations; however, some of the concepts from 
boundary-layer theory indicate the means to simplify the equations. 

The only assumption which can be made is the predominance of the convection in 
one main-flow direction. This leads to the (Reynolds number dependent) conclusion that 
diffusion can be neglected in this direction when compared with convection. Assiuning 
the x-direction to be the convective direction, the Navier-Stokes equations become 

for X -momentum: 



for y-momentum: 



for z -momentum: 



where the terms in braces are the streamwise x-direction diffusion ternis to be neglected; 
In addition, since the expression 

A = ^ + ^ + ^ ■ 

, ■ ■ ax 8y 8z 

also introduces x -diffusion terms in the y- and z -crossflow momentum eqiiations, these 
derivatives are also neglected here. The x-momentum equation is the same as would 
have been produced by a true boundary-layer scaling, but since no quantitative assump- 
tions have been made concerning the relative sizes of the x-, y-, or z-gradients, the 
y- and z-momentiun equations are retained, although in somewhat simpler form. 
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and the convective derivative ^ is given by 

■^ = U -^ + V -|- + w 
Dt 9x 3y 3z 

All x-gradients in the dissipation have also been neglected. Finally, the continuity equa- 
tion remains unchanged: . 

|;(PU) + ^V) + ^PW) = 0 . ( 5 ) 

^When these equations are supplemented by an equation of state and a viscosity relation-, 

P = pRT (6a) 


JL - 


where R is the imiversal gas constant, a system of five equations for five unknowns Is 
obtained after elimination of the density by the perfect gas equation of state. 

The elliptic nature of the Navier -Stokes equations In the x -direction has thus been 
eliminated; consequently, the equations are parabolic in x and marching integration 
may be used in the streamwise direction. This is significant computationally since it 
^ eliminates the need to store all the field quantities at each x-location which results in a 
substantial reduction in computer storage. Thus, the name parabolic -elliptic Navier- 
Stokes equations, since the assumptions allow a march in x away from kn-initial data 
plane, yet retain the elliptic character of the crossflow planes (Y-Z planes) due to the 
inclusion of all second derivatives in . y. and z. Flows with swirl or possible cross- 
flow recirculation (vortices) . in the Y-Z planes can be computed, and only reverse flow 
in the main -stream direction is eliminated due to the omission of x-diffusion. 

It should be noted that the continuity equation is hyperbolic; however, a march in 
the streamwise direction is still possible since the x -derivative can be expanded as 
follows: 
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= — p + 2- Uv 

T *^x X * 


Pu T 


and the term can be used to advance to the next station. There is no explicit diffu- 
sion term present and discontinuities which may be present can be expected to persist 
for large x-distances; hence, smooth initial profiles are desirable. 

Equations similar to these have recently been used for supersonic flow past a sharp 
cone at incidence (ref. 3) and for hypersonic leading -e^e flows, where a more formal 
analysis and scaling can be invoked (ref. 4). To date, no detailed work has been published 
on three-dimensional free -mixing flows. 


Validity for Supersonic Jet Flows 


Free jet and wake flows are common aerodynamic phenomena. These flows are 
generally turbulent, and the calculation of two-dimensional or axisymmetric turbulent 
free jets or wakes is difficult (ref. 2) because of problems associated with turbulence 
modeling; higher order modeling (two-equation models) is necessary in many cases. For 
three-dimensional flows with two essential cross-plane velocities, very few calculations 
have been made. To assess the modeling procedures for a three-dimensional flow, a cal- 
culation procedure valid for laminar flows, preferably in primitive variables to allow 
ease of incorporation of the turbulence models, is required. The parabolic -elliptic equa- 
tions (eqs. (1) to (6)) need to be verified for laminar jet calculations prior to their appli- 
cation to turbulent flows. 


If the equations are cast in nondimensional form by using free -stream values of 
Uoo, Pool and some suitable reference length which characterizes the problem, 

the equations become 


for X -momentum: 


nDu 
^ Dt 




yMc 


2 ax R Vy 


. a ,,/au\ 


(7) 


for y -momentum: 


^ Dt 




yMc 


2 ay R ]dy 



+ .2- 


ay 3 

az 

^[dzl 


( 8 ) 


for z -momentum: 


oSe = 
^ Dt 


1 9P 

yMoo^ 



( 9 ) 
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for energy: 



, for continuity: 


and for Sutherland's law: 

M = (^)t3/2 . (12) 

Rubin and Lin (ref. 5) have shown that equations of this t 3 ^e are singular at M = 1 
if the term is treated exactly and singular at M = 0 if the term is calculated 
in an explicit manner during the numerical calculation. U the p^ term is neglected 
entirely or specified from a boundary-layer approximation, then the parabolic march in 
X can proceed without difficulty. Thus, it was felt that since the p^ term should be 
included if necessary, the problem chosen to test the overall method should avoid any of 
these obvious integration difficulties. The M » 0 behavior takes place near boundaries 
(where u - 0), and therefore, the free jet problem avoids this singular behavior. How- 
ever, the jet caimot exhaust into an ambient atmosphere since here too u = 0. Thus a 
jet issuing into a moving free stream is suitable. To avoid any difficulties at M = 1 
both streams were chosen to be svq)ersonic. Thus, equations (7) to (12) will be solved 
for a supersonic jet issuing into a supersonic free stream. 

NUMERICAL PROCEDURE 
Integration Technique 

An implicit numerical procedure was chosen to solve the governing equations for a 
number of reasons. The success of implicit methods on the two- and three-dimensional 
boundary -layer equations implies that they should be efficient for the boimdaiy -layer - 
like parabolic -elliptic Navier -Stokes equations (eqs. (7) to (12)). It is e}q)ected that solu- 
tions will be required at large distances downstream from the initial data plane; conse- 
quently, large x-steps are desirable. The need to eliminate the step size restrictions of 
explicit methods leads to a consideration of xmconditionally stable methodis which are con- 
sistent in their marching variation. The particular implicit method used in this study is 
the alternating -direction Implicit (ADI) method of Peaceman and Rachford (ref. 6). The 
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ADI method is ideally suited for the solution of equations (7) to (12): There is no stabil* 
ity restriction on the step size, and hence, large x -steps are permitted. The meWc^ has 
second-order truncation error in its marching variation, which is also a requirement for 
the type of flow envisioned, since the x-history of the flow must be traced accurately at 
each step. Finally, the method does not require the inversion of a sparse -banded matrix, 
as a fully implicit or Crank -Nicolson method would.: Simple tridiagonal coefficient : 
matrices are generated at each step which require much less storage and time for their 
inversion in relation to sparse-banded matrices. The method has previously been shown 
to be effective for a set of equations similar to those used in the present approach (ref. 7). 

The ADI procedure is used to difference in x, with the y- and z -derivatives 
expressed as central differences, with the option of a nonunUorm grid included; that is, f 


« i _ au _ 

®2^,k = az - 


(Az+)(Az_)[(Az+) + (Az_)] 


,2ui ■ A, „ 

3z^ I (Az+)(Az_)[(Az+) + (Az.)] 


An example of the complete differencing scheme is shown for the x -momentum equation: 
as follows: 






f 

i+ir 

Ifl 

f£] ^ 1. 

ma' 



i+4 


, 4 . 


1+4 


U 2g 2„i . l/'M') 

'^j,k °y “j,k * 2\dT). 


j,k 


1+4 

% ^6-u.* t 

\9z/j,k 


1+4 


J,k 




i+i 

*“j.k«k2u,*k, 


(13a) 
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,i+l „♦ ^ 

*i,k -Pj,W 

Ax/2 j 


I r 1+^ l+'S 1 + i ' - . O • *'*■« 

4 . 1 }^(^^^\ ^ 2 i+i /Q^\ 2 1+1 25 2„i+l . if^lA ^ /aTV 2 

♦RWirjj.k |(arji.k °y“i.k * «yTj,tj .- t‘i,k h %k - 2(-5x)j_fe 


* ©,_; '»'^i*k "J,k “j*k 


(13b) 


The only difficulty arises from the cross derivatives of velocity present in the y- 
and z -momentum equations. These cannot be differenced implicitly since the tridiagonal 
structure of the resulting matrices would.be destroyed. These are treated in the same 
manner as all the nonlinear coefficients present in the differenced equations. 

Linearization Scheme 

The ADI procedure, which is second^order accurate, centers the x-derivatives about 
the point ^l + (See fig. 2.) This point is not equivalent to the Intermediate step of the 
ADI procedure. Hence, a method must be developed to compute all the nonlinear coeffi- 
cients (and the cross derivatives) at the ^l + point. This can be accomplished by a 

quadratic extrapolation from the two previous x-stations (i) and (i - 1) (ref. 7) or by use 
of the predictor -corrector procedure of Douglas and Jones (ref. 8 ). Both of these require 
additional storage. The procedure used in the present work handles the nonlinear terms 
by an iterative technique similar to that used for boundary -layer calculations in which a 
Crank -Nlcolson Integration is used (ref. 9). Any coefficient Q is calculated at the mid- 
point ^ simple formula 


,«i-k^4H:k*«i,k) 


The value at (I + 1) is not known on the first iteration and so the value is used. After 
the Integration to (I + 1) has been completed, new values of are computed, and the 

Integration step to (I + 1) is repeated to either a specified convergence or a specified 
number of iterations. At convergence this yields second-order -accurate coefficients 
which match the accuracy of the integration. The cross derivatives are also treated this 
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Solution of Matrix Equation 


The resulting linear difference equations are a set of five equations in five unknowns. 
Rather thanj make the quite arbitrary choice of the order of solution if a sequential teclT- ' 
nique of solving each equation in turn were elected, the five equations are solved simvil- • 
taneo\isly. The resulting matrix equation has a block tridiagonal structure which can be 
represented as 



where the unknown vector Wj contains all the five unknowns 

Wj = (uj,Vj,Wj,Tj,pj)'^ ( 16 ) 

the coefficients [a], [^, and [c] are simply the matrix coefficients of each particular 
imknown; that is, 



and the vector D is a source term in each equation. The components amhj represent 
the coefficient of the nth unknown from the mth equation at point j. 

The inversion of this block matrix is particularly simple. It consists of rewriting 
the usual tridiagonal algorithm with all multiplications replaced by matrix multiples, and 
all divisions replaced by matrix inversions. This simultaneous solution procedure was 
used previously by Krause, Hirschel, and Bothmann (ref. 9) where the pressure,- of course, 
was not one of the unknown variables due to the boundary-layer assumptions. ' The proce- 
dure is quite advantageous because, in addition to eliminating; the previously mentioned 
choice of solution order, it models the physics more precisely by allowing changes in any 
variable to be instantaneously sensed by all the others. It is also believed that this pro- 
cedure aids in the convergence of the iteration since it eliminates the use of lagged infor- 
mation previously calculated in a sequential solution. Thus it is fully implicit in the 
sense that a sequential solution is much like a Gauss -Seidel iteration which is explicit. 

The number of grid points used for the calculation was governed by conflicting 
requirements, ^or , acceptable resolution, many points were desirable in thes shear layer 
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between the jet and free stream and in the -jet itself. Also, reasonable distances away 
from the high shear region were necessai*y to allow boundary conditions in the free 
stream to be applied without distorting the results interior to the computational domain. 
However, excessive mesh points result in imacceptable machine storage requirements, 
increase computational time extravagantly, ^d impair job turnaroimd time. An initial 
compromise was to use a 41 x 41 grid to compute the quarter -plane of the jet flow, using 
the symmetry axes of the jet as boundaries. An equally spaced grid of Ay = Az = 0,1 
was first utilized, thereby placing a jet half -width at five grid spacings away from the 
axis and the outer edge of the domain seven times beyond this (40 grid spacings away 
from the axis). The required computer storage was 130 OOO 3 on the C DC 6600 system. 

Initial Conditions 

The initial conditions at the data plane representing the orifice location x = 0 
were chosen in a very rudimentary manner. Since too many points would be necessary 
to describe the merging of the free stream and duct flows just past a jet exit, and storage 
limitations were severe enough before this consideration, the initial velocity profiles 
were chosen as shown in figure 3. The jet and free stream are represented by two dis- 
tinct inviscid flows separated by a sharp boundary. Computationally this yields a one- 
grid-point discontinuity between uj^t and Uoo. This initial condition is probably the 
most severe that can be imposed while still generating an eventually realistic flow 
description. The initial conditions on the crossflow velocities were also modeled simply 
and were set equal to zero. 

The pressure distribution was chosen to be uniform at the free -stream level since 
an unmatched static pressure would undoubtedly produce shocks. These were not con- 
sciously sought as part of the problem, and the initial conditions were set to try to avoid 
their generation. The streamwise velocity and temperature levels were computed by 
assuming constant total temperature in the jet and free stream and specif 3 dng the jet and 
free -stream Mach numbers. 

. The free jet problem was expected to encounter few boundary condition troubles due 
to the avoidance of all u = 0 boundary conditions. The conditions placed on the variables 
were symmetry with no crossflow on the axes of-the jet and consistent free-stream con- 
ditions. Difficulties encoimtered with the application of these conditions are discussed in 
the next section. 

RESULTS AND DISCUSSION 

A standard test case was chosen to check the numerical procedure: the free- 
stream Mach number was set equal to 5.0; the jet Mach number at the initial station was 
7.5. The reference length was set equal to the minimum initial orifice 'width of a 
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1x1 square jet. . The Reynolds number based on this reference length was set equal to 
lo3, and the Prandtl number was set equal to 1. 


Effect of Initial Conditions • 

The initial behavior (small x) for all cases computed regardless of the grid config- 
uration or boundary conditions applied was essentially the same. The discontinuous 
velocity and temperature profiles were smoothed out by the diffusive terms, but the con- 
tinuity equation reacted to these discontinuities differently. The initial and subsequent 
u, T, and p profiles are shown qualitatively in figure 4. Referring to the figure, in 
the region below the initial discontinuity, Ux < 0 and Tx > 0, and above, ux > 0 arid 
Tx < 0. The X -gradients in the continuity equation appear as follows: 


P 

T 


3u 

ax 


pu 


— + (y 

2 ax 


and z 


gradients) = 


(18) 


Above the original discontinuity then, p^j < 0, and below, Pj^ > 0, neglecting the y- and 
z -gradients which are smaller here than the x -gradients. Hence, the pressure develops 
a blip around the initial discontinuity (see fig. 4) which persists for some distance before 
the profiles lose the influence of the initial conditions. This high-pressure-gradient pro- 
file causes divergence of the normal velocities about the initial breakpoint until the 
entrainment -induced bovmdary-layer-like velocities away from the jet axis are established. 
Eventually the pressure profiles smooth out to the expected constant case, but since the 
continuity equation is hyperbolic and contains no damping, ripples in the pressure profile 
of 0.2 percent of the free stream are commonplace. However, this initial behavior is of 
no concern, except in its influence on the downstream results, since the parabolic approx- 
imation is not valid m this regiori of very high x-gradients. 

Since the pressure is expected to be approximately constant in the developed jet, an 
attempt to artificially drive the pressure to its constant value more quickly, was made. 

An artificial diffusion term was introduced into the continuity equation in the hope of 
quickly smoothing the pressure profile. However, it was found that any introduction of 
these terms generated more diffusion in the other variables as well. Only values of an 
artificial viscosity larger than the actual flow viscosity had any significant effect and no 
improvement was effected by the incorporation of these fictitious terms. Thus, the flow 
was permitted to. naturally adjust to the initial discontinuity. If this was smoothed ini- 
tially over a few grid points, the pressure disturbance was smaller and shorter lived. 
Consequently it appears that no difficulties will be encountered if correct initial data are 
prescribed from experiments or boundary-layer calculations. 
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Problems Associated With Application of Boundary Conditions 


The e^qjectation that the elimination of solid walls and stagnant regions would ease 
boimdary condition difficulties was not realized. The boundaries of the computational 
domain were the axes of symmetry of the jet and the point in the free stream considered 
to approximate infinity. The conditions to be applied at these axes were symmetry of all 
quantities except the cross -plane velocity normal to the axis which was equal to zero. 
Thus at y = 0, for example, the conditions to be applied were Uy = Wy = Ty = py = 0 
and V = 0. Writing the governing equations differenced on the boundary (see fig. 5) and 
using the above relations to eliminate the unknown points outside the domain created by 
the central differences resulted in the values of the functions on the boundary becoming 
an additional row of unknowns to the vector Wj in the block tridiagonal system. The 
inversion technique thus allowed the simultaneous solution of the boundary values as well 
as the interior points, and the resulting solution was smooth for all variables except the 
pressiire. At the symmetry plane the pressure profile contained a spike (see fig. 6) and 
the other unknoums had gradients which differed from zero. This entire difficulty was 
eliminated by imposing a quadratic fit for the zero gradient condition onto each of the 
variables. The relation 


'^1 


^^"2 - ^3 
3 


(19) 


Insures a zero gradient with second-order accuracy, which matches the truncation of the 
difference scheme, and retains the tridiagonal aspect of the solution matrices. ' , 


The "infinity” boundary conditions must be imposed at boundaries in the free 
stream. These boundary conditions could be Imposed at a suitable distance from' the 
symmetry plane if asymptotic conditions could be derived far from the jet. However, 
asymptotic conditions are not known for the three-dimensional jet, and even if they were, 
the question arises as to a suitable distance at which they could be applied. Treatment 
of these boundary conditions can be classified into two groups: alterations in the actual 
conditions imposed at the last grid point and, concurrent with these, changes in the grid 
size. 

At first, all variables were specified at their free-stream values; that is, 
u = T = p = 1 and v = w = 0 at the last grid point in the domain with a uniform grid. 

The calcula,tlon proceeded smoothly, but after the usual jet crossflow velocity established 
itself, a sharp, one -grid-width gradient developed at the outer edge of the computational 
domain due to the difference between the negative entrainment velocity and the iniposed 
zero value. (See fig. 7.) , To alleviate this condition, an expanding nonuniform grid was 
used in place of the imiform one. Where the uniform 41 x 41 grid had y = z = 4 as the 
last grid point, the nonuniform grid using the same number of points allowed the edge to 
be displaced to y = z = 17.5, after remaining uniform to y = z = 1. This resulted.in a 
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delay in the appearance of this boundary gradient until the solution marched further 
downstream; however, it did not eliminate the problem. 

The e^anding jet both deflects the crossflow away from the jet axes and entrains^ 
fluid from the free stream. An interaction between the jet and variables such as v, w, 
and possibly p is quite likely to be present, and consequently, the specification of values 
for these- quantities at some finite near point would be of doubtful validity. Instead, the 
equations should determine the necessary boundary effects. To do this, a quadratic 
extrapolation consistent with the second-order differencing was imposed on v, w, and 
p by using both the uniform and nonimiform grids. The pressure solution quickly deterio- 
rated and the cause was traced to the lack of actual specification of a pressure level. 

With zero gradient conditions on the axes and extrapolation at the outer edge, no fixed 
pressure value is specified. Hence, the extrapolation was limited to the crossflow veloc- 
ities, which proved more successful, and the x-marching proceeded to a distance of 
approximately 20 jet widths. However, difficulties again occurred at the outer edge of 
the domain as shown in figure 8. The expected crossflow velocities appeared and the 
extrapolation did not affect the profiles. But, as the jet proceeds downstream, it-grows 
and the zero velocity point of the crossflow moves out from the axes. Eventually, as 
shown in figure 8, this point moved completely out of the computational domain, leaving 
only outflow from the centerline which is not characteristic of a jet. 

Until this point was reached, the development of the profiles showed the correct, 
trend for the streamwlse velocity (see fig. 9) but, as later calculations indicated, the 
greatest x-dlstance achieved was well ahead of the end of the jet core region, even if the 
nonuniform grid was used. For the square jet initial profile, all the variables were com- 
puted symmetrically about the jet axis, and the pressure became uniform to within 0.2 per- 
cent of the free-stream value. - ' ■ - 

All these difficulties with edge conditions arise due to the growth of the jet with x 
when it . is computed by using the unsealed, spatial coordinates. To overcome this prob- 
lem area a transformation was introduced in the Y-Z plane to contain the jet totally ' 
within the computational domain. 


Transformation of Coordinates 


.The objective in transforming the coordinates is to aid in the imposition of "infinity” 
bovmdary conditions. The uncertainty in these. conditions can be overcome only if the 
exact free-stream conditions can be imposed. To enable this to be accomplished, the 
point at infinity must be mapped to a finite point in the transformed space. One of the 

1 . 

simpler transforniations by which this can be accomplished is ' 



Ay 

1 + Ay 


(20a) 


556 



(20b) 


r = Bz 
^ 1 +fiz 

which maps zero onto zero and infinity in the physical plane onto unity in the 77 -^ plane. 
With equations (20) the governing equations are as follows (with I » x) : 

puu| +. [pv + 2Ap(l - 77 )]a (1 - 77 )^u^ + [pw + 2Bp(l - ?)]b( 1 - C)^^ = Pt 

+ i[a2(1 - n)^{iir^T^Urj + u^) + b2(1 - ?)2(,x.j.t^u^ + u^^) (21) 

puv^ + jpv + 1 Ap(l - t?)1a(1 - 7j)2v^ + |pw + 2Bp(l - 0 ]b(1 - = - - 2 

•L _ J yMoo 

+ IJ-YjjjjY+ b2(1 _ C)^(MtT^v^ + 

+ AB(1 - 77)2(1 - ^)2 1 - I T^w^\ \ (22) 


puw^ + [pv + 2Ap(l - T7)]a(1 - ?7)2wj, + |pw + 1 Bp(l - C)1b(1 - = ^ B(1 - ^)2p^ 

' L • J ' . y *^00 

+ + I b 2(1 - ^)4^7i,j.T^w^ + pw^^) 


+ AB(1 - 77)2(1 - ^)2 1 + Pt(t77V^ - | 


puT^ + pv + ^ A(1 - 77)j A(1 - v)^Trj + pw + ^ B(1 - C)^ B(1 - ?)2 t^ = ^[up^ 

+ vA(l - 77)2p^ + wB(l - C)2p^j + a2(1; - 77)^[ptT^ 2 +,pT^^j + b2(1 - ^)4[ji,j,t^2 

+ pT^^)j + pi | a2(i . ^)4(u^2 + w^2 + 4 ^^2) b2(i . ^)4|^i^2 + v^2 + 4 ^^2j 

+ AB(1 - 77)2(1 . 5)2^w^v^ - I VjjW^j (24) 
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UP| - puT^ +.pu^ + AU - - pvTjj + pv^j + B(1 - C)^(wp^ - pwT^ + pw^j = 0 (25) 

These equations were cast in finite difference form and solved numerically by using the 
same ADI method, linearization scheme, and block tridiagonal solution algorithm as pre- ’ 
viously discussed. The boundary conditions u = T = p = 1 and v = w = 0 are imposed 
at 77 = 1 or C = 1» corresponding to infinity in the untransformed plane. A coarse 
21 X 21 equally spaced mesh was used oh the unit square in the 77 -^ plane to reduce the 
core storajge to more acceptable levels (eSOOOg) for computation. 

Results Computed in Transformed Plane ' 

Cases comparable to the ones calculated on the untransformed variables were com- 
puted successfully for a variety of initial geometries with no difficulty. The loss of 
entrainment velocities never occurred, confirming the usefulness of the transformation. 

The first test case was a square jet with unit sides. The streamwise centerline 
velocity decay of this jet is presented in figure 10. After a lengthy core region the cen- 
terline velocity quickly decays according to the laminar relation for axisymmetric jets 
x"l (ref. 1). Thus, although the equations compute the flow as if it, were truly three- 
dimensional, the axisymmetry is reproduced. Another square jet was computed to test 
the transformation. A large jet with sides of length equal to four was formed by taking . 
unequal scalings of y and z. That is, five points in C . sufficed to give z = 2 (half 
the jet width), while 11 points in 77 were necessary to give the equivalent y = 2. The 
centerline decay for this case is also shown in figure 10. After a much longer core than 
the 1 X 1 jet (since the initial shear layer is four times farther from the axis), this con- 
figuration also very rapidly begins to decay along the axisymmetric curve. In fact, if the 
1 X 1 curve is displaced to the right by the difference of core lengths, the decay curves 
coincide. 

The tluree-dimensional capabilities of the method were tested on rectangular jets of 
varying aspect ratios. A rectangular 2 x 1 jet gave results for velocity decay shown in 
figure 10. The initial core region was of the same length as the 1 x 1 square jet. ‘ 
However, the decay was slower than the lx 1 jet, and in fact there was a region where 
the two-dimensional laminar jet decay x"^/^ described the flow. Eventually though, 
the decay increased and approached the axisymmetric behavior seen before with slope 
x"l. Only this far -field behavior could be described by an axisymmetric boundary -layer 
analysis, whereas the length to its asymptotic decay could not. 

A rectangular 4 x 1 jet is more two-dimensional than either of the previous cases, 
as indicated by its decay. (See fig. 10.) ^ After the initial core length, characteristic of 
the distance needed for disturbances one width away from the centerline to reach the axis, 
the decay curve obviously follows the two-dimensional slope for a greater distance down- 
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stream to about x = 500 (see fig. 10) where it also gradually starts to approach the 
axisymmetric curve. The length x = 500 is not too surprising, for this is the core dis- 
tance of the 4x4 square jet, the distance needed for disturbances four widths away to 
reach the axis. The even higher aspect ratio 8 x 1 jet continued the established trend 
and is also shown in figure 10. 

; Finally,.- note the advantage gained through the use of the transformation. Prior to ; 
its use, by usir^ a greatly expanding mesh in the untransformed plane with extrapolation 
at the last mesh point, complete loss of entrainment took place at approximately x = 20. 
Even at this point the accuracy of the solution is questionable, and this is not even half- 
way to the end of the core region as computed by using the transformation and depicted 
in figure 10. 

CONCLUDING REMARKS 

The parabolic -elliptic Navier -Stokes equations have been shown to be a viable ' 
method for computation of three-dimensional supersonic jet flows. The difficulties asso- 
ciated with the unbounded domain of the jet can be eliminated by incorporating a transfor m 
mation into the equations so that the points of infinity in the cross plane are mapped^to a 
finite value at the transformed plane. 

Although the character of the computed flow can be predominantly axisi^metric in ' 
the case of the square jet or approach a two-dimensional limit, as is the case for the 
8x1 rectangular jet, no prior assumptions to this effect are required. The solution 
computes the particular flow vinder investigation as though it were three-dimensional, 
allowing the initial geometry of the prescribed jet to determine the ultimate nature of 
the solution. 

The present calculations have been for laminar jets; however, jet flows are gener- 
ally turbulent. Thus, even if some simpler method could be used instead of the present 
integration scheme to recover the gross characteristic of the laminar flow, the inclusion 
of turbulence modeling in the fully elliptic crossplane requires numerical treatment. 

A subsonic analog of the parabolic -elliptic Navier -Stokes equations, where the 
streamwise pressure gradient is correctly accounted for, is currently being investigated 
so that the more advanced state of the higher order incompressible turbulence models can 
be drawn upon for inclusion into the governing equations. If this proves successful, then 
the s\Q)ersonic eqimtions will be used as a means to test various modeling procedures for 
compressible turbulent flows. 
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(a) Upper surface blowing. (b) F-15 with proposed rectangular nozzle. 

Figure 1.- Examples of three-dimensional jets in aerodynamic flows. 
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Figure 2.- The ADI procedure: first step implicit in z, explicit in y; second step 


implicit in y, explicit in z. 
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Figure 5. - Differencing used on a boundary. 



attributed to boundary differencing. 
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Figure 7.- Crossflow profiles when 
free-stream conditions imposed: 



Figure 8. - Crossflow profiles when free-stream conditions extrapolated. 
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Flgfure 9.- Computed velocity profiles at various stations. 
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